\documentclass[border=10pt,tikz]{standalone}
\usepackage{pgfplots,pgfplotstable,tikz}
\usetikzlibrary{intersections,matrix}
\usepackage{verbatim}
\usepackage{moresize}
\usetikzlibrary{calc}
\usepackage{datatool}
\usetikzlibrary{shapes}
\usetikzlibrary{plotmarks}
\usepackage{csvsimple}
\usepackage{siunitx}
\usepackage{mathtools}
\usepackage{amsmath,amssymb}
\usetikzlibrary{shapes.geometric}
\usepackage{booktabs}
\usepackage{xspace} 
\usepackage{etoolbox}
\usepackage{multirow}
\usetikzlibrary{arrows.meta}
\usetikzlibrary{positioning}
\usetikzlibrary{decorations.text}
\usetikzlibrary{decorations.pathmorphing}
\pgfplotsset{compat=1.13}
\usepackage{fp}
\usepgfplotslibrary{fillbetween}

\definecolor{lightblue}{RGB}{76, 100, 237}
\definecolor{blue}{RGB}{19,28,84}


% Defines "sgn" function
\def\sgn{\mathop{\rm sgn}\nolimits} 
% Defines the expectation function
\def\E{\mathop{\mathrm{E}}\nolimits} 

% Defines a function that ignores the zero on the y- or the x-axis, such that the origin
% is only labeled once. Activate in the axis options by adding the key: "ignore zero=x" or ""ignore zero=y."
\pgfplotsset{ignore zero/.style={%
		#1ticklabel={\ifdim\tick pt=0pt \else\pgfmathprintnumber{\tick}\fi}
	}}
	
	
	
	\begin{document}
		
	% TAKE INPUTS FROM THE FOLLOWING FILE:
	% "diff_in_diff_ohie.xls" --- generated in Stata by the program file: "ohie_sumstats.do" 	(for plotting the "OR line of best fit")
	\pgfmathsetmacro\pB{0.15}
	\pgfmathsetmacro\pI{0.41}
	\pgfmathsetmacro\xcom{\pB*0.5+\pI*0.5} %midpoint of OR compliers
	\pgfmathsetmacro\xnt{0.5+\pI*0.5}	%midpoint of OR never takers
	\pgfmathsetmacro\SRHcomOR{0.546102642684361}	%value of "mean_UC" in the "selfhealth_fair" workbook of the excel spreadsheet
	\pgfmathsetmacro\SRHntOR{0.343792633015007}	%value of "mean_NT" in the "selfhealth_fair" workbook of the excel spreadsheet
	\pgfmathsetmacro\slopeOR{(\SRHntOR-\SRHcomOR)/(\xnt - \xcom)} %slope of the OR line of best fit 
	\pgfmathsetmacro\interceptOR{\SRHcomOR - \slopeOR*\xcom} %intercept of the OR line of best fit 
	
	% "diff_in_diff_brfss.xls" --- generated in Stata by the program file: "brfss_diff_diff.do" (for plotting the "MA line of best fit")
	\pgfmathsetmacro\pBMA{0.90}
	\pgfmathsetmacro\pIMA{0.95}
	\pgfmathsetmacro\xcomMA{\pBMA*0.5+\pIMA*0.5}	%midpoint of MA compliers
	\pgfmathsetmacro\xntMA{0.5+\pIMA*0.5}	%midpoint of MA never takers
	\pgfmathsetmacro\SRHcomMA{0.20594047254236}	%value of "mean_UC" in the "selfhealth_fair" workbook of the excel spreadsheet
	\pgfmathsetmacro\SRHntMA{0.175694292886576}	%value of "mean_NT" in the "selfhealth_fair" workbook of the excel spreadsheet
	\pgfmathsetmacro\slopeMA{(\SRHntMA-\SRHcomMA)/(\xntMA - \xcomMA)} %slope of the MA line of best fit 
	\pgfmathsetmacro\interceptMA{\SRHcomMA - \slopeMA*\xcomMA} %intercept of the MA line of best fit 
	
	% [yaxismin, yaxismax, yaxislength] Calculate graphing region for the y-axis
			\pgfmathsetmacro\yaxismin{0}
			\pgfmathsetmacro\yaxismax{0.82}
			\pgfmathsetmacro\yaxislength{\yaxismax-\yaxismin}
			
		
			
			
			
			% MTE(pC) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
			\begin{tikzpicture}[>=stealth,element/.style={align=left},square/.style={regular polygon,regular polygon sides=4}]
			% Axis options
			\begin{axis}[					
			width=15cm,					% width of chart
			height=10cm,				% height of chart
			xlabel={$U_D$:  unobserved net cost of treatment, Oregon},				% x-axis title
			ylabel={Fair or Poor Health, Untreated},					% y-axis title
			ylabel style={rotate=0},  % rotate the y-axis title
			samples=100,  				% sets the resolution for curved lines
			xmin=0,xmax=1,			    % x-axis range
			ymin=\yaxismin,
			%ymin=1.2*\Dinterceptr,ymax=1.05*\yaxismax,			% y-axis range
			ignore zero=y,			    % ignore the zero on the y-axis (defined in preamble)
			xtick={0,\pB,\pI,1},	    % only show these x ticks
			xticklabels={0,$p_{C}^{OR}=0.15$, $p_{I}^{OR}=0.41$,1}, % label x ticks
			ytick={0, \SRHcomOR, \SRHcomMA, \SRHntOR, \SRHntMA},	    % only show these y ticks
			ymax = \yaxismax,
			legend style={draw=none, at={(0.5,\yaxismax+0.12)},anchor=north},	% set legend style, specifically no outline
			legend cell align={left},	% align legend content to the lfet
			clip=false,					% place nodes outside of plot area
			xticklabel style={
						align=center,
						yshift=-5.5ex,
						name=label,
						fill=none,
						append after command={(label.north) edge [solid] ++(0,6ex)},
										/pgf/number format/fixed,
										/pgf/number format/precision=2},
					%x label style={at={(axis description cs:0.5,-0.2)},anchor=north}, % ncrease distance between title and x-axis
						yticklabel style={	/pgf/number format/fixed,
						/pgf/number format/precision=2},
					% place some extra y ticks ...
					extra y ticks={0},
					% ... but do not repeat the already existent labels
					% (from the "normal" ticks) ...
					extra y tick labels={},
					% ... and set the style of the extra ticks to `major'
					extra tick style={
						grid=major,
					},
			%ticklabel style={fill=white},
			%xlabel style={align=center,xshift=0pt,yshift=-20pt},
			 mydot/.style={
			 	circle,
			 	fill=white,
			 	draw,
			 	outer sep=0pt,
			 	inner sep=1.5pt
			 },
			]		
			
			%draw extra minor ticks manually
			\draw [thin,draw=gray] (\pIMA,\yaxismin) -- (\pIMA,\yaxismin+0.014);
			\draw [thin,draw=gray] (\pBMA,\yaxismin) -- (\pBMA,\yaxismin+0.014);
			\draw [thin,draw=gray] (\pIMA,\yaxismax) -- (\pIMA,\yaxismax-0.014);
			\draw [thin,draw=gray] (\pBMA,\yaxismax) -- (\pBMA,\yaxismax-0.014);

			% Dashed Vertical Lines OR
			\addplot[
				line width=1pt, black, dashed, opacity=100, thin
			] coordinates {(\pB,0) (\pB,\interceptOR + \slopeOR*\pB)};
			
			\addplot[
				line width=1pt, black, dashed, opacity=100, thin
			] coordinates {(\pI,0) (\pI,\interceptOR + \slopeOR*\pI)};
			
			
			%Dashed Vertical Lines MA
			\addplot[
				line width=1pt, black, dashed, opacity=100, thin
			] coordinates {(\pBMA,\interceptMA + \slopeMA*\pBMA) (\pBMA,\yaxismax)};
	
			\addplot[
				line width=1pt, black, dashed, opacity=100, thin
			] coordinates {(\pIMA,\interceptMA + \slopeMA*\pIMA) (\pIMA,\yaxismax)};
			
												
			% OR SRH line
			\addplot[
				line width=1.8pt, blue, dashed, dash pattern=on 3.8pt off 2pt
			] coordinates {(0,\interceptOR) (1,\interceptOR + \slopeOR)}{};
			
			\node[outer sep=1.3pt,square, fill=blue,inner sep=2pt] (P) at (\xcom, \SRHcomOR) {};
			\node[outer sep=1.3pt,square, fill=blue,inner sep=2pt] (P) at (\xnt, \SRHntOR) {};
		
		
			% MA SRH line
			\addplot[
				line width=1.8pt, lightblue, dashed, dash pattern=on 3.8pt off 2pt
			] coordinates {(0,\interceptMA) (1,\interceptMA + \slopeMA)}{};
				\node[outer sep=1.3pt,square,fill=lightblue,inner sep=2pt] (P) at (\xcomMA, \SRHcomMA) {};
			\node[outer sep=1.3pt,square,fill=lightblue,inner sep=2pt] (P) at (\xntMA, \SRHntMA) {};
						
			
			\addplot[opacity=0, name path=Z] coordinates {(0,0) (1,0)};
			
	
						
				\node[label={[align=center, fill=none]-270:{\small $p_{I}^{MA}$ \\ $\pIMA$}},circle,fill=none,inner sep=2pt,opacity=0] at (axis cs:\pIMA+0.005,\yaxismax+0.067) {};
						\draw [thin,black] (\pIMA,\yaxismax) -- (\pIMA,\yaxismax+0.067);
						
						\node[label={[align=center,fill=none]-270:{\small $p_{C}^{MA}$ \\ $\pBMA$}},circle,fill=none,inner sep=2pt,opacity=0] at (axis cs:\pBMA-0.005,\yaxismax+0.067) {};
						\draw [thin,black] (\pBMA,\yaxismax) -- (\pBMA,\yaxismax+0.067);
						
						\node[label={[align=left,fill=none]-270:{0}},circle,fill,inner sep=2pt,opacity=0] at (axis cs:0,\yaxismax+0.067) {};
						\draw [thin,black] (0,\yaxismax) -- (0,\yaxismax+0.067);
						
						\node[label={[align=left,fill=none]-270:{1}},circle,fill,inner sep=2pt,opacity=0] at (axis cs:1,\yaxismax+0.067) {};
						\draw [thin,black] (1,\yaxismax) -- (1,\yaxismax+0.067);
						

						
			
			% add legend

			\legend{,,,,{\small OR line of best fit},{\small MA line of best fit}}
			
			
		
			%Put two Complier group labels on x-axis
			\draw[draw=none,decorate,decoration={brace,amplitude=10pt,mirror,raise=4pt},yshift=0pt] (\pB,\yaxismin) -- (\pI, \yaxismin) node[pos=0.5, below, yshift = 0pt] {Compliers};		
			
			\draw[draw=none,decorate,decoration={brace,amplitude=10pt,mirror,raise=4pt},yshift=0pt] (0,0) -- (\pB, 0) node[pos=0.5, below, yshift = 0pt, align=center] {Always \\ Takers};
			
			\draw[draw=none,decorate,decoration={brace,amplitude=10pt,mirror,raise=4pt},yshift=0pt] (\pI,0) -- (1,0) node[pos=0.5, below, yshift = 0pt] {Never Takers};
			
\draw[draw=none,decorate,decoration={brace,amplitude=10pt,raise=4pt},yshift=0pt] (0,\yaxismax) -- (\pBMA,\yaxismax) node[pos=0.5, above, yshift = 0pt] {Always Takers};

\draw[draw=none,decorate,decoration={brace,amplitude=10pt,raise=4pt},yshift=0pt] (\pBMA,\yaxismax) -- (\pIMA,\yaxismax) node[pos=0.5, above, yshift = 0pt] {C};

\draw[draw=none,decorate,decoration={brace,amplitude=10pt,raise=4pt},yshift=0pt] (\pIMA,\yaxismax) -- (1,\yaxismax) node[pos=0.5, above, yshift = 0pt] {NT};
		
		
			\node[label={[align=center,fill=none]-270:{$U_D$: unobserved net cost of treatment, Massachusetts}},circle,fill,inner sep=2pt,opacity=0] at (axis cs:0.5,\yaxismax+0.11) {};
				

			\end{axis}
			


			
			\end{tikzpicture}
			
		
	
		
	\end{document}
